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Abstract 

Nonhnear feedbacks in the Earth System provide mechanisms that 
can prove very useful in understanding complex dynamics with rela- 
tively simple concepts. For example, the temperature and the ice cover 
of the planet are linked in a positive feedback which gives birth to mul- 
tiple equilibria for some values of the solar constant: fully ice-covered 
Earth, ice-free Earth and an intermediate unstable solution. In this 
study, we show an analogy between a classical dynamical system ap- 
proach to this problem and a Maximum Entropy Production (MEP) 
principle view, and we suggest a glimpse on how to reconcile MEP 
with the time evolution of a variable. It enables us in particular to 
resolve the question of the stability of the entropy production maxima. 
We also compare the surface heat flux obtained with MEP and with 
the bulk-aerodynamic formula. 

1 Introduction 

A very broad class of problems in climate modelling consists of studying the 
evolution of a particular field (e.g. surface temperature, precipitation, etc) 
when an external parameter, or forcing^ is varied. Most of the time, the 
response to this variation is not linear. Feedbacks can amplify or damp the 
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effect of the initial perturbation. One of these feedbacks aroused a proficient 
branch in scientific hterature in the 70s', when Budyko and Sellers simultane- 
ously suggested that the interaction between sea ice and climate could have 
dramatic consequences. Indeed, the higher the global temperature on Earth, 
the less the ice cover is likely to extend, and thus the lower the albedo. A 
lower albedo leads in turn to a higher global temperature, and so on and so 
forth until all the ice is melted. Stimulated by this pioneer work, the ques- 
tions of the stability of the climate as well as the consequences such feedbacks 
might have for understanding paleoclimates were extensively studied, using 
the whole hierarchy of models, from the most simple Energy Balance Models 
(EBMs) to the complex General Circulation Models (GCMs). 

Using ID EBMs, Budyko and Sellers had found two stable equilibrium 
positions for the edge of the ice cover, one corresponding to the present 
climate and one to a fully ice-covered Earth ^ i55j. A large part of the 
subsequent work was concerned with verifying that these results still held 
with various different versions of the Budyko-Sellers models, with different 
heat transport par ameterizat ions, temperature dependance expressions in 
the planetary albedo, numerical schemes,... [TOl Ell ESI [39l e.g.]. Some 
elegant analytical solutions were found for these models O |39l [40j , and var- 
ious mathematical methods were applied to determine the stability of the 
equilibria [191 EHl ESI SI US] . Owing to the extreme sensitivity of climate to 
variations in the solar constant found by the first studies, the precise posi- 
tion of the tipping point between present climate and a deep freezed Earth 
was of primary concern. Further investigation by [29] and [43 J revealed that 
the sensitivity was much less than initially thought. A fundamental ques- 
tion raised by these results was that of the transitivity of the climate system 
in Lorenz's terminology [301 [3T] . and the difference between forced and free 
fiuctuations [54l [191 El- For a comprehensive review of the various models, 
parameterizations and problems pertaining to Energy Balance Models and 
the ice-albedo feedback, the reader is referred to [41j. 

In this contribution, we will first give a brief account of the reformulation 
of these questions with the vocabulary of dynamical system theory: how do 
multiple equilibria arise from the ice-albedo feedback, what does the bifurca- 
tion diagram look like, etc. The model used here is a two box energy balance 
model with a simplified radiative transfer using the Net Exchange Formu- 
lation (see e.g. [9j), and a bulk aerodynamic formula for the surface heat 
fiux. In a second step, we draw an analogy between this dynamical system 
view and the results obtained when predicting the surface heat fiux from 
the Maximum Entropy Production (MEP) principle. The MEP principle, as 
originally expressed by [46] [471 HH] for the climate system, provides a varia- 
tional principle to compute energy fiuxes that are not otherwise constrained 
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by the laws of physics. OriginaUy, Paltridge and others apphed MEP to the 
meridional energy transport [46l [471 EOj [HJ |32l ^-g-]? but other studies [441I52] 
indicate that it may be valid on the vertical also. 

As noticed by [43], [6] and [T4j, the bifurcation giving birth to multiple 
equilibria in the case of the ice-albedo feedback has a fundamentally radiative 
nature, and has nothing to do with transport properties of the atmosphere. 
This encourages one in thinking that a zero-dimensionnal model is sufficient 
to capture the structure of the mechanism while avoiding the use of more 
cumbersome mathematics (namely the Sturm-Liouville theory, required for 
one-dimensional models such as [19j). Therefore we will restrict ourselves 
here to this idealized case. Note also that most of our work could be trans- 
posed easily to other feedbacks, like the water- vapour feedback. 

2 The ice-albedo feedback, multiple equilibria 
and the hysteresis cycle: the dynamical sys- 
tem approach 

2.1 A simple two-layer EBM using the Net Exchange 
Formulation 

We use a slightly different formulation of the model described in [24j, as 
represented in Fig. [TJ A grid cell is characterized by a surface temperature 
Tg and an atmospheric temperature T^, and we note (respectively ^^Y) 
the flux of solar energy received by the ground (respectively absorbed by 
the atmosphere). Radiative exchange use the Net Exchange Formulation^ 
in which the basic objects are not energy fluxes at a given level but rather 
the energy exchange rate between two layers in the atmosphere or between 
one layer and a boundary surface (see [9]). is the net energy exchange 
rate between the ground and the atmospheric column per unit surface (ie the 
greenhouse effect), and (respectively ^If) is the cooling to space term 
for the atmosphere (respectively the surface). The net energy exchange rates 
per unit surface are expressed as functions of Tg and as: 



*r = {s{c^)-s){l-a)iS, (1) 

vl>ff = {s + as*)iS, (2) 

K = taT^-t^T^^ (3) 

"^i^ = taT^. (4) 

K - (l - ^) (5) 
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where a is the Stefan-Boltzmann constant, a is the surface albedo, t, 5, 5*, 5 
are radiative coefficients, S is the solar constant, ^ accounts for the annual 
mean zenith angle of the sun and ji is the Elsasser factor (see [24j for a 
derivation of the equations and a discussion of the coefficients). 

In addition to radiation, energy is exchanged due to atmospheric and 
oceanic transport as well as surface heat ffuxes. Let us merge all these en- 
ergy transfer modes into two variables: 7^ (respectively 7^) represents the 
net convergence (the opposite of the divergence) of energy into the atmo- 
spheric cell (respectively the surface layer). Writing C^a for the atmospheric 
convergence (this variable was designated by C in |24j), Co for the oceanic 
convergence (this was not taken into account in [21]), and q for the surface 
heat ffux, we have 

la = Ca + q^ (6) 

Ig = Co-q- (7) 

Knowing the convergence of energy in each cell - atmosphere or ground - 
it is in general not possible without further assumptions to separate the con- 
tribution due to surface ffuxes, atmospheric transport, and oceanic transport 
when applicable. Of course, over land, it is reasonable to assume that 7^ is 
just the surface energy ffux (ie = 0) , and then 7^ + 7^ is the convergence 
of energy due to the atmospheric ffow. In this study, as we will only use the 
zero-dimensional version of this model, we will always have (a = Co = 0, and 
thus 7« = -7^ = q. 

At steady-state, the energy balance equations for the atmosphere and the 
surface read 

i?a(T„T,)+7, = 0, (8) 

i?,(r„r,) + 7, = 0, (9) 



where 



IR 



(10) 



= (s + as*)^S + t{aT;-2aT^) 
and 

^^g\-^a) -^g) ^ gs ^ ag ^ sg 

= {s -s){l- a)iS - t{aT^ - aT^) - (l - ^) ^^^^ 

In this form, the steady-state equations cannot be solved since 7^ 

and 7^ are unkown. In the next section we introduce a parameterization of 
these quantities as functions of and Tg. In section [3| we use the MEP 
principle to compute them. 
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2.2 The zero-dimensional model with bulk aerodynamic 
formula 



In the case of a zero-dimensional, two-layer model considered here, the net 
convergence of energy in the atmospheric box (ie the divergence of the dia- 
batic heating at the surface, = Q = ~lg) can be simply interpreted as the 
surface heat flux. In this section, we adopt a bulk aerodynamic formula 
to express this flux as a function of the temperatures and Tg\ 

la = qbaf{Ta, Tg) = CpaCoUs (T, - T,) . (12) 

where Cd '^^ the drag coeflicient, Ug is the surface wind speed and c^a is the 
heat capacity of the atmosphere per unit surface area (similarly c^g is the heat 
capacity of the ground). Now the model can be seen as a two-dimensional 
dynamical system: 



with 



T 

^9 



F(T„T,), (13) 



F(T„.r,)=( 1. (14) 



and 



Fi(r„,r,) = —{RaiTa,T,) + q,afiTa,Tg)), (15) 



F2{Ta,T,) = ±{R^(Ta,T,)-q,af{Ta,Tg)). (16) 



Cpg 



Our main interest here is to flnd the equilibrium positions of the system, 
ie the flxed points of the dynamical system, given by the roots of F, and to 
study their stability. Of course, the dynamics of a two-dimensional dynamical 
system can be more complex than just a relaxation to an equilibrium position 
(although it is still rather gentle, see [22] for example), contrary to one- 
dimensional dynamical systems. Still, let us note here that the flrst equation 
in F{Ta^Tg) = can be solved algebraically in to obtain a relation = 
f{Tg) where (T^^Tg) is a flxed point of the system. Thus the number of 
flxed points of the two-dimensional system is exactly the number of roots of 
the scalar equation F2(/(T^),T^) = 0. 

For simplicity, we will consider here the projection of the dynamical sys- 



tem (13) onto the Tg axis: 

f, = F2ifiT,),T,). (17) 
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As just explained, this dynamical system, although not mathematically 



equivalent to the full dynamical system (14), has the same equilibrium posi- 
tions. Physically, this simplification is motivated by the fact that the atmo- 
sphere can be assumed to reach equilibrium very quickly, hence the evolution 



of Ta is enslaved by the dynamics of Tg. In other words, the system (17) is 



just the system (14) with Cpa = 0. 



2.3 Multiple equilibria 

The values of the coefficients used here are reproduced in table [T] Taking for 
the albedo the fixed value ao = 0.15, the system only has one fixed point, 
as plotting the function F2{f{Tg)^Tg) = clearly shows (see Fig. [2]). In this 
case, the equilibrium is at a global mean surface temperature of 288 K. 

But in reality, the higher the global mean temperature, the lower the 
extent of the regions that can sustain an ice-cover. This positive feedback 
can be encoded in the following temperature dependance for the albedo : 

where (respectively aj) represents the value of the planetary albedo 
over an ice-free (respectively fully ice-covered) area, and Tq and AT are 
parameters determining the transition from ice-free to ice-covered conditions 
(see Fig. |3]). One could simply use a step function between ice- free and 
ice-covered albedo values, or a piecewise linear function, but we choose this 
expression because it depends smoothly on the temperature. 



Replacing a in Eq. (14) with (18) yields a new dynamical system 



T 

f 
^9 



G{Ta,Tg), (19) 



where the fixed points are again determined by the conditions, g being 
defined similarly to / (or obtained by substitution of the albedo function into 



T: = g{T;), (20) 

= G'2(^(r;),r;). (21) 



Plotting the curve G2{g{Tg)^ Tg) as a function of Tg (Fig. |4]) shows that for 
certain values of the solar constant, three solutions coexist. This range can 
be determined to be approximately 0.98aSo < S < I.OSaSq. Outside this range. 
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only one solution subsists. For the present value of the solar constant, S = Sq^ 
for instance, these equilibria correspond to a fully glaciated Earth {snowball 
state) ^ 249 K, an ice- free Earth ^ 287 K which can be identified 
with the present climate, and an intermediate glacial state Tg ^ 275K. For 
a low value of the solar constant (e.g. 0.955^0), only the snowball state Tg 
subsists. Similarly, at high solar constant (e.g. l.lS'o), the only equilibrium 
is found on the ice-free branch . 

A fixed point X* of the dynamical system X = F[X) is said to be (lin- 
early) stable if all the eigenvalues of the jacobian of F are negative (see [Ij for 
a complete classification of the two-dimensional fixed points). In this model 
we find that Tg and Tg are always stable nodes when they exist, while T^ 
is a saddle-point. 

The stability can also be read directly on Fig. [4] for the ID-reduced 
system: stable equilibria correspond to roots of the function with negative 
derivative, while at the unstable equilibrium, the curve crosses the x-axis 
with an upward slope. 

Summarizing the above results. Fig. [5] represents the curve of the fixed 
points when sweeping a large range of values for S : it is the bifurcation 
diagram of the dynamical system. Creation of a pair of stable/unstable 
equilibria at the tipping points 0.985^0 and l.OS/So is called a saddle-node 
bifurcation. Thus the hysteresis curve obtained for the temperature stems 
from the bifurcation structure of the dynamical system as two back-to-back 
saddle-node bifurcations. It is noteworthy that this figure does not depend 
upon the particular coefficients choice in the bulk formula, nor on the green- 
house effect. Would we set g^^/ = (radiative equilibrium with greenhouse 
effect) or/and t = (greenhouse effect shut down), the hysteresis curve would 
remain qualitatively the same. 



2.4 Potential for the dynamical system 



The full two-variables dynamical system ( |13| ) cannot be expressed as the 
gradient of a potential function, but its one-dimensional projection can, like 
any other one-dimensional dynamical system. Let us thus introduce the 
potential V (defined up to an additive constant) such that 

dV 

Fixed points of this dynamical system correspond to critical points (ie 
extrema in this ID case) of the potential. The stability criterion becomes 
that stable fixed points are minima of the potential: 
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- ^ < 0' (23) 
while its maxima are unstable fixed points. 

Figure [6] shows the shape of the potential for different values of the solar 
constant. At low solar constant (e.g. 0.955^0), the potential has only one 
critical point, a minimum at T ?^ 245X. Increasing the value of the solar 
constant levels down the potential curve, until a second local minimum ap- 
pears (along with a local maximum) with T above the freezing point, around 
S ^ 0.985^0. At = S'o, it is clear that the potential has two minima at 
T ^ 250K and T ^ 290K and a maximum at T 275K. Further increase 
of the solar constant leads to a deeper minimum at T > ° C while the 
minimum at T < ° C becomes shallower. Around S ^ I.OS/Sq, the mini- 
mum at T < ° C disappears (it annihilates with the local maximum) ; for 
S = l.lS'o, the only minimum is found at T 300 K. 

Note that, as expected, the critical points of the potential obtained for 
the three values of the solar constant considered here match with the values 
of Fig. [4j Also, the number of critical points of the potential changes at the 
bifurcation points of the dynamical system. 



3 The entropy production rate and the ice-albedo 
feedback 

In this section, we do not use anymore the bulk aerodynamic formula for the 
surface flux 7^ = —7^, but the Maximum Entropy Production Principle, as 
described in [24j. The first application of the MEP principle to climate is 
found in [46], where the meridional energy transport in a zonally averaged 
box-model is chosen so as to maximize the entropy production. The resulting 
climate is in striking accordance with observations. In spite of successful ap- 
plications in other areas as well, the domain of validity of the MEP principle 
remains unclear due to the lack of a fully convincing proof (see [71 18] and the 
comments in [211 [2]). More details on theoretical issues and practical use can 
be found in [45l EH |35] . 

3.1 The entropy production rate in zero-dimensions 



Let us consider the model of Sect. |2.1| and introduce the entropy production 
rate per unit surface 

^ = ^ + ^- (24) 
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Substituting Eqs. g-g into Eq. Q for 7„ and 7^, (J can be considered 
as a functional of the temperature field. We are looking for its maxima 
subject to the constraint 

la+lg = 0. (25) 

The sum of equations ([s]) and ^ can thus be solved for Ta as a function 
of Tg^ and the entropy production rate a is simply a function of one variable. 
Its graphical representation for the set of parameters given in table [l] (fixed 
albedo ao) is shown in Fig. [Tj It is clear that there is only one local maximum, 
corresponding to a surface temperature Tg ^ 295K. 

Now, replacing in the equations the constant albedo ao by the temperature- 



dependent albedo (18), the resulting entropy production rate curve is plotted 



in Fig. [8] for different values of the solar constant. 

Unlike the potential for the dynamical system in Sect. 2^, the entropy 
production rate always has at least two local maxima and a local minimum. 
In fact, over a rather narrow range, estimated to be 0.95aSo < S < 1.005aSo, 
the entropy production rate even has three maxima and two minima. This is 
even clearer on the contour plot of the entropy production rate as a function 
of Tg and S/Sq (Fig. [9]). Hence, there is indeed an analogue of the fold of 
the potential in the classical dynamical system picture in the context of the 
entropy production surface, but the values at which it takes place do not 
exactly correspond. 

Besides, a large portion of the curve on Fig. [8] lies under the abscissa 
axis: for the corresponding range of temperature values, the entropy produc- 
tion rate is negative, contrary to what the second law of thermodynamics 
states (or more precisely its extension to non-equilibrium systems). It seems 
reasonable to impose the condition 

a(T,) > 0, (26) 

thereby restricting the range of values Tg can actually take. In this case, 
this is equivalent to requiring that the surface heat flux goes from hot to 
cold. With this additional constraint, the range of possible values of the 
solar constant allowing for coexistence of multiple equilibria (two or three) 
can be determined approximately: 0.8aSo < S < 1.125'o. 



3.2 Stability of the MEP states 

In the classical understanding of the MEP principle, the rate of entropy pro- 
duction a is a function deflned on the manifold of steady-states which reaches 
a maximum at the most probable state. In the presence of several local max- 
ima, it is generally believed that the final equilibrium state of the system 
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will be the global maximum. In our case, there are three local maxima of 
the entropy production rate for the present value of the solar constant, as 
discussed in the previous section. We know from the dynamical system ap- 
proach that there can indeed be several steady-states (that coincide with the 
positions of the entropy production maxima, as discussed in the previous sec- 
tion) for a given set of parameters, and the actual steady-state of the system 
is determined from the initial conditions. In the absence of fluctuations, the 
system remains in this state. Hence it is certainly not suflicient to retain 
the global maximum of the entropy production rate as representing the flnal 
state of the system. Instead one must flnd a practical way to select a local 
maximum for given initial conditions. As a particular case, we would obtain 
a way to distinguish between local entropy production maxima representing 
dynamically stable steady-states and dynamically unstable ones. 

This involves the introduction of time in the MEP formulation. So far, 
there was no mention of time in the MEP approach as we were only concerned 
with steady-states. Even though we claim that the entropy production max- 
ima correspond to equilibrium points, —a is by no means a potential for the 
dynamical system. Indeed, the dynamics of the system is simply given by 
the first law of thermodynamics. Here, it reads 

Cpa^ = RaiTa,Tg)+^a, (27) 

dT 

cp.^ = RATa,T,)+jg. (28) 

Similarly to the steady-state entropy production rate, we can define the 
instantaneous entropy production rate: 

' ~ Ta{t) ^ T,{t) ~Tait) dt 



(29) 



using Eqs. (27)-(28). Note that the instantaneous entropy production rate 
ai and the steady-state entropy production rate a coincide at steady-state. 

As ai appears as the natural generalization of a taking into account the 
time derivative of the dynamical variables, we suggest that the system may 
follow the trajectory maximizing the instantaneous entropy production rate, 
seen as a function of the time-dependent unknown fiuxes 7^, jg (always linked 
by the relation ja + 7g = 0). This approach is very similar to what [25] 
advocates for. 

In practice, it is easier to reformulate the above suggestion with a time- 



discretized system (see Fig. 10). Let us consider two snapshots of the system 
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separated by a finite time interval dt. We note T^,Tg the values of the air 
and surface temperature at time t. The instantaneous entropy production 
rate becomes: 

— - RaiTlTl)^ + - [c,, ' - R,{TlTl)j 

(30) 

Suppose we know the state of the system at time t — 1 (ie ^ and Tg ^ are 
given) . Then our postulate is that and can be chosen so as to maximize 
aj (with fixed T^~^ and Tg~^) subject to the constraint 7^ + 7^ = 0. Iterat- 
ing this process leads to a trajectory maximizing the instantaneous entropy 
production rate at each timestep, starting from a given initial condition. 

Integrating the system with this method, initialized in the vicinity of the 
different maxima of the entropy production rate at steady-state, provides a 
criterion for stability: it is found here that the warm branch as well as the 
snowball branch of Fig. |9]are stable, while the intermediate branch is unsta- 
ble. The maxima of the entropy production and their stability are plotted as 
functions of the solar constant on Fig. [11} analogously to Fig. [5| This result 
draws the final line in the parallel between the dynamical system approach 
and the MEP approach. Note that the limits of this analogy are reached at 
some points: figure [TT] cannot be considered as a usual bifurcation diagram. 
As a consequence, the lines of existence of the maxima need not depend 
continuously on the parameter, and for certain values of the parameter (for 
example S ^ O.OS'o), two stable maxima coexist with no unstable manifold 
to separate them. 

The trajectory maximizing the instantaneous entropy production rate 
in the way explained above thus yields stability properties for the differ- 
ent steady-states that are consistent with the dynamical system approach. 
Hence, it seems legitimate to use this hypothesis as a relaxation equation^ 
in a similar fashion as [53j. However, there is no certainty that the system 
actually follows this maximum instantaneous entropy production trajectory. 
It would be very valuable to investigate the range of validity of this new 
application of the MEP principle in future studies, theoretically or on other 
examples. We can already adduce some material to support our relaxation 
equations approach. In fact, the only novelty as compared to the common 
use of MEP in the steady-state context is the inclusion of time derivatives of 
the dynamical variables in the entropy production rate. But one can simply 
consider these time derivatives as known fiuxes, playing exactly the same 
role as Ra{Ta^ Tg) or Rg{Ta^ Tg). The only difference is that computing these 
fiuxes requires that we consider a bigger system (here simply the state of the 
system at times t — 1 and t), even though the number of unknowns in the 
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big system remains the same (T^ and Tj, whereas T^~^ and Tg~^ are fixed). 
In this respect, there is no fundamental difference between the time dimen- 
sion and any geometric dimension, which are customarily included in MEP 
models. 

Alternatively, one could consider the total entropy production rate (ie 
the integral of the instantaneous entropy production rate over time) as a 
functional of trajectories and claim that the system follows the trajectory that 
maximizes this functional subject to the relevant constraints ([l2l[T3j, [llj 
and ^36j have developed this idea in the case of Markov chains by maximizing 
the information entropy as a function of both the probability of the states 
and that of the transition rates). As we should show in a forthcoming study, 
this is particularly suitable for periodic phenomena, such as the seasonal 
cycle. Regarding the stability of the steady-states, we expect this method 
to yield the same results as the maximum instantaneous entropy production 
relaxation used here. 



3.3 Surface heat flux and SnowbaU Earth deglaciation 

In the case of the first section, the surface heat fiux is parameterized as a 
function of and Tg. As a consequence of this strong constraint, one could 
draw a bifurcation diagram for g^^/ very similar to Fig. [5j with relatively 
weak surface heat fiux for low solar constants (around 20 W.m~'^)^ strong 
surface heat fiux g^^j at high solar constants (around lOOVI^.m"^), with an 
unstable branch linking the two. 

On the contrary, the surface heat fiux obtained through the MEP proce- 



dure Qmep is much less constrained by the temperature gradient. Figure 12 
shows the surface heat fiux as a function of the temperature gradient Tg — Ta 
for both cases: q^af and qmep- It is clear that the two differ completely, not 
only because the temperature gradients in the different climates are very dif- 
ferent, but also because the shape of Qj^ep as a function of the temperature 
gradient is far from linear. Note that in the MEP snowball state, although 
the temperature gradient is relatively high, the surface fiux remains very low. 
On the warm branch for the MEP state, high values of qj^ep are obtained for 
high values of the solar constant. Hence, decreasing the solar constant brings 
the surface fiux down, until the point where only the snowball state survives, 
with a similar low value of the surface heat fiux. 

This discrepancy between the two graphs is likely to be significant: it has 
been suggested that the suppression of the vertical temperature gradient in 
the snowball state numbers amongst the reasons that make deglaciation of 
the snowball Earth so difiicult [50[[5ll[23. Indeed, the temperature inversion 
isolates the surface from all the forms of energy exchange: the greenhouse 
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effect can only warm the surface when the air aloft is colder, latent heat 
plays a very limited role in this very dry atmosphere, and the sensible heat 
flux is also restricted by the vertical structure of the atmosphere. [50j points 
out that a crucial role may be played by the surface fluxes parameterization 
and the convection parameterization. Here the simplicity of the model does 
not allow us to discuss the static stability, nor to come up with a clear 
explanation of the questioning Fig. [12} but it does certainly reinforce the idea 
that surface heat fluxes parameterization can play critical parts on important 
paleoclimate problems. In the case of the MEP surface heat flux, our results 
tend to indicate that it would be possible for the snowball earth to withstand 
a vertical temperature gradient higher than expected with very little loss in 
the form of sensible heat, thereby damaging the thermal shield of the surface 
layer mentioned above. 

On a similar note, f34] performed a thorough investigation of the thermo- 
dynamic properties of the snowball Earth as compared to warm climates in 
the model of intermediate complexity PLASIM [15], using the formalism of 
non-equilibrium thermodynamics applied to the climate system as described 
in [33j. Computation of the thermodynamics efliciency, irreversibility and 
material entropy production clearly characterizes distinct thermodynamic 
regimes for the snowball Earth and ice-free climate. Our remarks about 
the surface heat flux in snowball conditions add up to their thermodynamic 
analysis. 



4 Conclusion 

The analogy developed in this study leads to some enlightening conclusions. 
First, about the ice-albedo feedback in itself, it provides a variational princi- 
ple different from those previously suggested, with a thermodynamic motiva- 
tion. On the contrary, all the candidates for variational formulations of the 
problem examined previously were rather ad hoc potentials for the dynam- 
ical system. The parallel between potentials properly speaking, which fully 
describe the dynamics of the system, and the entropy production rate, which 
only characterize equilibrium states, was pushed one step further with the 
introduction of a method to integrate a trajectory using the MEP principle. 
In particular we have shown that this method predicts the correct stability 
for the MEP predicted equilibria. We also investigated the behaviour of the 
surface heat flux in the snowball state. The results hint that MEP might 
prove useful in such extreme situations where the usual parameterizations 
face important difliculties. However, the highly simplifled model considered 
here does not allow us to conclude against or in favour of the MEP parame- 
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terization, as compared to the bulk- aero dynamic formula. 

As far as the MEP conjecture is concerned, our work adds up to the rel- 
atively short list of efforts up to now (essentially |56l [57] and [26]) to sort 
out how the principle should be understood in the presence of multiple en- 
tropy production maxima. [57j suggested that a dynamical system, in their 
case the thermohaline circulation, when multiple steady-states are available, 
should move to the most dissipative one. [371 EH] showed strong limitations 
to this interpretation in full generality. Here, we find that steady-states of 
a system with unknown turbulent fiuxes correspond to local maxima of the 
entropy production seen as a function of the unknown fiuxes. The stability 
of these maxima does not seem to depend on the numeric value of the en- 
tropy production at that point. Instead, we suggest that the question of the 
dynamic stability can be investigated by a relaxation process maximizing the 
instantaneous entropy production rate. 
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Table 1: Values for the parameters of the OD model (radiative coefficients, 
bulk aerodynamic formula parameters and ice-albedo feedback parameteri- 
zation) . Note that the values for the heat capacities depend on the thickness 
of the layer and on the nature of the surface (ocean or land) , but this has no 
influence on the steady-state results. 
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Figure 1: A grid cell of the model, adapted from [24]. "i/^j are the energy 
exchange rates per unit surface due to radiative transfer (see text), q is the 
surface heat flux and (a is the atmospheric energy convergence. Over the 
oceans, there is also an oceanic energy convergence (o- 
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Figure 2: Function F2{f{Tg),Tg) as a function of Tg (see text) with a fixed 
albedo has only one root. 
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Figure 3: Surface albedo a as a function of surface temperature Tg in K. 
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Figure 4: Function G2{g{Tg),Tg) as a function of Tg (see text) including the 
ice-albedo feedback for different values of the solar constant: 0. 955*0 (blue), 
(red), l.OSS'o (yellow), 1.15*0 (green). 
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Figure 5: Bifurcation diagram of the bulk aerodynamic formula model. Tq^ 
Tp and Ts are plotted against S/ S{) when they exist. Stable fixed points are 
plotted in blue while the unstable solution is in dotted red. This figure clearly 
shows that two saddle-node bifurcations occur at respectively S ^ 0.985^0 and 
S ^ 1.08^0. 
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Figure 6: Potential V (normalized) as a function of temperature Tg (in K) 
for three different values of the solar constant: 0. 955^0 (red), Sq (blue) and 
1.1 S'o (yellow). For the present value of the solar constant, the potential has 
a double well shape, with two stable equilibria, while for the two other values, 
the potential has only one minimum. 
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Figure 7: Entropy production rate as a function of the surface temperature 
Tg for the CD model ai S = Sq. The only local maximum corresponds to 
Tg ^ 295K. 
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Figure 8: Entropy production rate as a function of the surface temperature 
Tg for the OD model with ice-albedo feedback. For a low value of the solar 
constant {S = O.SS'o, blue curve), there is only one local maximum with 
positive entropy production rate. The same holds for high solar constant {S = 
1.25^0, yellow curve), while there are three local maxima and two minima, all 
with positive entropy production rates, for S = Sq (red curve). 
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Figure 9: Contour plot of the entropy production rate as a function of the 
solar constant (normalized by its present-day value) and the surface temper- 
ature Tg (in K). Negative contour lines are dashed, positive contour lines are 
solid and the null contour line is the thick solid line. Shades of blue represent 
negative values of the entropy production rate. 
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Figure 10: To discuss the stability of the steady states predicted by MEP, we 
need to extend the principle to obtain a time-dependent formulation. This is 
done by maximizing the instantaneous entropy production rate. To compute 
the time derivative of the temperature, we consider it as a known flux in time 
seen as a geometric dimension of the space upon which MEP operates (see 
text). In green, the fluxes that can be computed from the state variables 
(T^^-^T^^Tj-^Tj). in red, the unknown flux obeying MEP. 
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Figure 11: Entropy Production maxima as a function of the solar constant, 
normalized by its present value. The solid lines (respectively the dotted line) 
correspond to dynamically stable (respectively unstable) equilibria in the 
sense of paragraph |3.2[ Note that this is not a bifurcation diagram in the 
usual meaning. 
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Figure 12: Comparison between the bulk aerodynamic formula surface heat 
flux (top) and the MEP predicted surface heat flux (bottom) as a function of 
the temperature gradient Tg — Ta. The red solid line corresponds to the warm 
branch of the bifurcation diagram, the blue solid line to the snowball state 
and the dotted yellow line is the unstable branch. Note the very different 
scales for — T^. 
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